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(19.  continued) 

- - ->>  Rayleigh  waves  normally  incident  upon  2-D  shallow  heterogeneity  are  simulated  by  the  linear  finite- 

difference  method  to  study  attenuation,  transmission,  and  reflection  of  Rayleigh  waves  and  to  measure  the 
Rayleigh-to-P  and  -SV  body  wave  conversion.  Transmission,  reflection,  and  scattering  depend  on  the  depth, 
average  scale  size  of  the  heterogeneity,  and  the  amplitude  of  the  spatial  fluctuation  of  velocity.  As  expected, 
larger  spatial  variation  in  velocity  attenuates  Rayleigh  waves  more  than  smooth  media,  and  the  attenuation  is 
roughly  proportional  to  the  variance  of  the  velocity  fluctuation.  The  attenuation  and  scattering  due  to  shallow 
heterogeneity  is  weaker  than  attenuation  due  to  moderately  rough  topography. 

Scattered  body  wave  energy  is  studied  as  a  function  of  frequency,  scattering  angle,  and  wave  type  (P 
or  SV).  Attenuation  of  Rayleigh  waves  by  scattering  from  2-0  shallow  velocity  heterogeneity  is  dominated 
by  conversion  to  body  waves  and  in  particular  SV  energy.  Low  frequency  P  and  SV  energy  is  scattered  in  a 
backwards  direction,  and  high  frequency  P  and  SV  energy  is  scattered  in  a  forward  direction. 

As  with  scattering  from  rough  topography,  much  of  the  converted  SV  energy  will  be  trapped  in  the 
crustal  waveguide  at  Lg  phase  velocities.  Therefore,  Rayleigh  (Rg)  to  SV  body  wave  conversion  by  shallow 
heterogeneity  and  topography  should  contribute  to  the  formation  of  Lg  by  explosions,  quany  blasts,  and  shal¬ 
low  earthquakes. 

A  comparison  was  made  with  results  for  P-coda  from  Greenfield  (1971).  The  comparison  indicates 
that  self-similar  and  Gaussian  models  could  be  derived  with  rms  velocity  variations  between  7  and  15%  in 
the  upper  3  km  of  the  crust  that  would  produce  the  observed  P-coda/P  power  levels  observed  by  Greenfield 
(1971). 

Linear  finite-difference  (FD)  method  was  used  to  compare  the  excitation  of  far-field  P-  and  SV-waves 
generated  by  shallow  dilatational  sources  in  a  suite  of  heterogeneous  2-D  crustal  models.  The  crustal  models 
tested  included  simple  layered  structures,  media  with  random  velocity  perturbations  having  Gaussian  or  self- 
similar  autocorrelation  functions,  media  with  rough  or  gentle  topography  generated  by  Markov  chains,  and 
laminated  media  with  sinusoidal  folds.  The  numerical  experiments  were  conducted  by  directing  a  broadband 
planar  P-  or  SV-wave  with  appropriate  incidence  angle  upon  the  testing  models.  Jhf  dilatational  strain  Jus-  - 
lory  at  a  shallow  linear  array  of  grid  points  was  then  recorded  so  that  the  far-field  P-  or  SV(Lg)-waves  from 
shallow  dilatational  sources  could  be  inferred  by  use  of  the  principle  of  reciprocity.  The  raw  FD  synthetics 
were  deconvolved  so  as  to  represent  the  response  due  to  explosion  sources  with  a  fixed  yield.  The  mean  peak 
amplitude  of  the  synthetics  for  each  model  are  compared  to  that  for  a  reference  model  consisting  of  a  simple 
layered  medium.  The  average  energy  content  in  an  appropriate  signal  window  was  measured  as  a  comple¬ 
ment  to  the  amplitude  measurement  Both  approaches  show  essentially  the  same  pattern  of  P/SV  excitation, 
namely  that  models  with  topography  consistently  produce  the  strongest  P-SV  conversion  among  all  types  of 
crustal  models.  The  introduction  of  interfaces  (e.g.,  dipping  layers)  alone  does  not  by  itself  increase  SV 
excitation  with  tire  required  slowness  range.  Thus  «*(P)  -  i%(Lg)  appears  to  be  smaller  for  models  with  topo¬ 
graphic  relief  (e.g.,  the  Degelen  region  of  the  central  portion  of  the  East  Kazakh  Test  Site  (EKTS))  than  for 
models  with  dipping  layers  or  folded  sedimentary  rocks  (e.g.,  Shag  an  River,  eastern  EKTS).  This  result  is 
quite  different  from  Nuttli’s  (1987)  observations  based  on  WWSSN  film  chip  readings  of  Lg.  which  suggest 
that  m*(P)  -  m^Lg)  varies  from  0.036  ±  0.01S  for  the  Shagan  River  area  to  0.27  1 0.03  for  the  Degelen  area. 

Recommendations  for  further  work  include: 

(1)  Extensions  of  the  current  finite-difference  code  from  2-D  to  3-D  to  study  the  attenuation  of  body 
waves  by  3-D  heterogeneity  in  the  crust,  test  hypotheses  about  the  generation  of  P  coda  and  anisotropic 
P  wave  generation,  and  generation  of  transverse  Lg  by  explosions. 

(2)  Introduction  of  other  numerical  methods  to  explore  the  coupling  (scattering)  of  modes  of  wave-guide 
regional  phases  such  as  Pg  and  Lg,  as  well  as  the  scattering  of  Pn  and  Sn.  These  methods  include  2-D 
and  3-D  scattering  from  localized  heterogeneity  as  well  as  from  rough  boundaries. 

(3)  Coupling  of  efficient  reflectivity  methods  to  finite  difference  calculations  to  propagate  the  scattered 
field  to  regional  distances  and  to  drive  the  finite  difference  responses  with  realistic  in-coming  regional 
phases. 

(4)  Investigation  of  scattering  of  fundamental  and  higher  mode  short-period  Rayleigh  waves  by  2-D  topog¬ 
raphy  and  shallow  heterogeneity  with  more  realistic  velocity  gradients  near  the  surface. 

(5)  Extension  of  the  general  topographic  boundary  condition  to  include  the  general  fluid-solid  interface  for 
the  modeling  of  scattering  at  rough  fluid-solid  boundaries. 

(6)  Improvement  of  the  polygonal  free-surface  boundary  conditions  for  higher  precision. 
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SUMMARY  OF  RESEARCH  COMPLETED  DURING  FEB  86  -  FEB  88 

FD  Simulations  of  Rayleigh  Wave  Scattering  by  Shallow  Heterogeneity 

Rayleigh  waves  normally  incident  upon  2-D  shallow  heterogeneity  are  simulated  by 
the  linear  finite-difference  method  to  study  attenuation,  transmission,  and  reflection  of 
Rayleigh  waves  and  to  measure  the  Rayleigh-to-P  and  -SV  body  wave  conversion  ( cf 
AFGL-TR-87-0322,  also  Figures  1,  2,  and  3).  Transmission,  reflection,  and  scattering 
depend  on  the  depth,  average  scale  size  of  the  heterogeneity  and  the  amplitude  of  the  spa¬ 
tial  fluctuation  of  velocity.  As  expected,  larger  spatial  variation  in  velocity  attenuates 
Rayleigh  waves  more  than  smooth  media,  and  the  attenuation  is  roughly  proportional  to 
the  variance  of  the  velocity  fluctuation  (Figures  4  and  S).  The  attenuation  and  scattering 
due  to  shallow  heterogeneity  is  weaker  than  attenuation  due  to  moderately  rough  topogra¬ 
phy. 

Scattered  body  wave  energy  is  studied  as  a  function  of  frequency,  scattering  angle, 
and  wave  type  (P  or  SV).  Attenuation  of  Rayleigh  waves  by  scattering  from  2-D  shallow 
velocity  heterogeneity  is  dominated  by  conversion  to  body  waves  and  in  particular  SV 
energy.  Low  frequency  P  and  SV  energy  is  scattered  in  a  backwards  direction,  and  high 
frequency  P  and  SV  energy  is  scattered  in  a  forward  direction. 

As  with  scattering  from  rough  topography,  much  of  the  converted  SV  energy  will  be 
trapped  in  the  crustal  waveguide  at  Lg  phase  velocities.  Therefore,  Rayleigh  (Rg)  to  SV 
body  wave  conversion  by  shallow  heterogeneity  and  topography  should  contribute  to  the 
formation  of  Lg  by  explosions,  quarry  blasts,  and  shallow  earthquakes. 
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A  comparison  is  made  with  results  for  P-coda  from  Greenfield  (1971).1  The  com¬ 
parison  indicates  that  self-similar  and  Gaussian  models  could  be  derived  with  rms  velo¬ 
city  variations  between  7  and  15%  in  the  upper  3  km  of  the  crust  that  would  produce  the 
observed  P-coda/P  power  levels  observed  by  Greenfield  (1971)  (Figure  6). 


1  Greenfield,  R.  J.  (1971),  Short-period  P-wave  generation  by  Rayleigh-wave  scattering  at  Noway*  Zeralya,  J.  Gtopkys.  Rtt.,  76, 
7988-8002. 
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Figure  4.  Attenuation  factor  observed  as  a  function  of  frequency  from  the  FD  simulations,  "x",  "+", 
triangles,  and  circles  correspond  to  fluctuation  of  P  wave  velocity  v  ■=  20%,  10%,  7%,  and  5%  respec¬ 
tively.  The  shallow  heterogeneous  media  tested  are  (from  top  to  bottom):  4  Gaussian  media  with  a  - 
lkm,  h  ■  3.2km;  4  Gaussian  media  with  a  -  2km,  h  *  3.2km;  4  self-similar  media  with  a  -  1km,  h  - 
3.2km;  4  self-similar  media  with  a  -  2km,  h  -  3.2km;  folded  sinusoidal  layers  with  h  ■  3.2km,  X  - 
2km,  peak-to-peak  amplitude  2.5km. 
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Figure  5.  The  attenuation  factor  1/Q  at  0.78  Hz  and  I.56Hz  versus  energy-flux  £  of  Gaussian  models 
of  various  thickness  (lkm,  2km,  and  3.2km)  and  velocity  fluctuations  (v  -  5%,  7%,  10%  and  20%).  (A) 
0.78Hz,  a-1,  fitted  to  curve  T  =  r^11®42,  0.78Hz,  a-2,  fitted  to  curve  T  =  r<£12991,  (C)  1.56  Hz,  a=l, 
fitted  to  curve  T  =  r.50M34,  (D)  1.56  Hz,  a-2,  fitted  to  curve  T  =  r^o  m44. 
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Figure  6.  Power  spectral  ratios  of  the  scattered  P  and  S  waves  to  the  incident  Rayleigh  wave  of  various 
models.  Units  are  3.4*  10-3  and  I.2*1(T3  erglseclcir?l{cm  of  heterogeneity)/(cm2  incident  Rayleigh  wave) 
for  the  P-wave  and  S-wave  coda  power  density  at  a  depth  of  3.2  km  in  the  grid.  P-wave  coda  is  lower  set 
of  values,  S-wave  coda  is  upper  set  of  values.  (A)  Gaussian  autocorrelation  models  (a-=  1  km,  h  -  3.2  km), 
v  varies  from  20%  (top),  10%,  7%,  to  3%  (bottom).  (B)  Same  as  (a)  except  a-  2  km.  (C)  Self-similar 
autocorrelation  models  with  a-  1  km.  (D)  Self-similar  autocorrelation  models  with  a-  2  km.  (E)  Folded 
models. 
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SUMMARY  OF  RESEARCH  COMPLETED  DURING  THE  PERIOD  MAY  87  TO  FEB  88 

FD  Studies  of  P-SY(Lg)  Coupling  in  2D  Crustal  Models 

A  linear  finite-difference  (FD)  method  has  been  used  to  compare  the  excitation  of  far- 
field  P-  and  SV -waves  generated  by  shallow  dilatational  sources  in  a  suite  of  hypothetical 
heterogeneous  2-D  crustal  models  (cf.  AFGL-TR-88-0025).  The  crustal  models  tested 
included  simple  layered  structures,  media  with  random  velocity  perturbations  having  Gaussian 
or  self-similar  autocorrelation  functions,  media  with  rough  or  gentle  topography  generated  by 
Markov  chains,  and  laminated  media  with  sinusoidal  folds.  The  numerical  experiments  were 
conducted  by  directing  a  planar  P-  or  SV-wave  with  appropriate  incidence  angle  upon  the  test¬ 
ing  models  (Figures  7  and  8).  The  dilatational  strain  history  at  a  shallow  linear  array  of  grid 
points  was  then  recorded  so  that  the  far-held  P-  and  SV(Lg)-waves  from  shallow  dilatational 
sources  could  be  inferred  by  using  the  principle  of  reciprocity.  The  raw  FD  synthetics  were 
deconvolved  so  as  to  represent  the  response  due  to  explosion  sources  with  a  fixed  yield  (Figure 
9).  The  mean  peak  amplitudes  of  the  synthetics  for  each  model  are  compared  to  that  for  a 
reference  model  consisting  of  a  simple  layered  medium.  The  average  energy  content  in  an 
appropriate  signal  window  was  measured  as  a  complement  to  the  amplitude  measurement  (Fig¬ 
ure  10).  Both  the  amplitude  and  energy  measurements  show  essentially  the  same  pattern  of 
P/SV  excitation,  namely  that  models  with  topography  consistently  produce  the  strongest  P-SV 
conversion  among  all  types  of  crustal  models  (Tables  1  and  2).  The  introduction  of  interfaces 
{e.g.,  dipping  layers)  alone  does  not  by  itself  increase  SV  excitation  with  the  required  slowness 
range.  Thus  mb(P)  -  mb(Lg)  appears  to  be  smaller  for  models  with  topographic  relief  than  for 
models  with  dipping  layers  or  folded  sedimentary  rocks. 

These  synthetic  results  are  consistent  with  observations  for  Novaya  Zemlya  (Nuttli, 
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1988)2  and  Shagan  River  (Nuttli,  1986), 3  baaed  on  WWSSN  film  chip  readings  of  Lg.  Novaya 
Zemlya,  which  has  rough  topography,  shows  relatively  higher  Lg  with  respect  to  P  (mb(P)  - 
nifcCLg)  -  -0.11)  than  does  the  somewhat  flatter  Shagan  River  test  site  (mb(P)  -  m^Lg)  -  0.04). 
However,  Nuttli  (1987)4  also  obtained  an  even  lower  value  of  Lg  relative  to  P  (mb(P)  -  m^Lg) 
-  0.27)  for  the  Degelen  Mountain  test  site,  only  70  km  away  from  Shagan  River.  If  this 
Degelen-Shagan  bias  is  teal,  then  it  must  be  due  to  near-source  effects,  and  these  cannot  be 
explained  by  the  FD  results  obtained  to  date.  However,  some  recently  archived  high-quality 
digital  seismograms  recorded  at  the  Chinese  Digital  Seismic  Network  indicate  more  Lg  excita¬ 
tion  (with  respect  to  P)  at  Degelen  than  at  Shagan  River  (Figure  13),  which  is  consistent  with 
the  numerical  results  (Figures  11  and  12). 

The  finite-difference  results  also  show  negatively  correlated  P  and  SV  energy  (Figures  1 1 
and  12),  which  provides  a  preliminary  explanation  of  the  success  of  the  unified  yield  estimator. 
Measuring  all  possible  phases  tends  to  reduce  the  effects  of  uneven  energy  release  on  source 
size  estimation.  To  understand  this  issue  in  a  more  quantitative  manner,  and  to  derive  an 
optimal  weighting  scheme  to  combine  all  phases,  theoretical  studies  with  numerical  simulations 
on  detailed  deterministic  (rather  than  oversimplified  or  hypothetical)  models  of  the  Soviet  test 
sites  are  necessary. 


1  Nuttli,  O.  W.  (1988),  Lg  magnitude*  and  yield  ettimate*  (or  underground  Novaya  Zemlya  auclear  explosions,  Bull.  Seam 
Soc.  Am.,  78,  873-884. 

5  Nuttli,  O.  W.  (1986),  Lg  magnitude*  of  selected  East  Kazakhstan  underground  explosions,  Bull.  Seam  Soc  Am.,  76, 
1241-1251. 

Nuttli,  O.  W.  (1987),  Lg  magnitudes  of  Degelen,  East  Kazakhstan,  underground  explosions,  Bull  Seam  See.  Am.,  77, 
679-681. 
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TABLE  1.  AMPLITUDE  COMPARISON  OF  P  AND  SV  Lg  EXCITATION 

model 

P 

Lg 

W!3 

Description  of  the  model 

0 

0.000 

0.000 

0.000 

reference  model  (1 -uniform  layer,  5+0%,  2km  thick) 

1 

-0.207 

0.202 

-0.409 

rough  TOPO  +  1  uniform  layer  (5+0%,  2km  thick) 

2 

-0.006 

0.132 

-0.139 

gentle  TOPO  +  self-similar  layer  (5+10%,  2km) 

3 

-0.196 

0.110 

-0.305 

rough  TOPO  +  Gaussian  layer  (5+10%,  2km) 

4 

-0.023 

0.073 

-0.096 

gentle  TOPO  +  1  uniform  layer  (5+0%,  2km) 

5 

-0.034 

0.044 

-0.078 

self-similar  layer  (5+10%,  2km  thick) 

6 

-0.162 

0.019 

-0.181 

folded  sinusoidal  layers  (L=2,H=2.5,5+20%) 

7 

-0.031 

0.014 

-0.045 

folded  sinusoidal  layers  (L=2,H=2.5,5+10%) 

8 

-0.134 

-0.037 

-0.098 

self-similar  layer  (5+20%,  2km  thick) 

9 

-0.029 

-0.037 

0.008 

folded  sinusoidal  layers  (L=5,H=2.5,5+10%) 

10 

0.003 

-0.058 

0.061 

2 -Gaussian  layer  (4.5+10%/5+10%,  total  2km) 

11 

0.019 

-0.091 

0.110 

steeply  dipping  layers  (52°) 

12 

0.011 

-0.093 

0.104 

gently  dipping  layers  (26°) 

13 

0.018 

-0.137 

0.155 

steeply  dipping  layers  (-52°) 

14 

0.009 

-0.143 

0.152 

gently  dipping  layers  (-26°) 

TABLE  2.  COMPARISON  OF  P  AND  SV  Lg  SPECTRAL  CONTENT  ON  0.5- 1.0  Hz 

model 

P 

Lg 

E5 S3 

Description  of  the  model 

■an 

0.000 

0.000 

0.000 

reference  model  (1 -uniform  layer,  5+0%,  2km  thick) 

1 

-0.400 

0.083 

-0.483 

rough  TOPO  +  uniform  layer  (5+0%,2km) 

2 

-0.049 

0.057 

-0.106 

gentle  TOPO  +  self-similar  layer(5+10%,2km) 

3 

-0.363 

0.063 

-0.426 

rough  TOPO  +  Gaussian  layer  (5.0+10%, 2km) 

4 

-0.180 

0.019 

-0.199 

gentle  TOPO  +  uniform  layer  (5+0%, 2km) 

5 

0.016 

0.009 

0.007 

self-similar  layer  (5+10%,2km) 

6 

0.099 

-0.031 

0.130 

folded  sinusoidal  layers  (5+20%,L-2,H«2.5) 

7 

0.058 

-0.101 

0.159 

folded  sinusoidal  layers  (5+10%,L=2,H=2.5) 

8 

-0.026 

-0.049 

0.023 

self-similar  layer  (5+20%, 2km) 

9 

0.015 

-0.163 

0.178 

folded  sinusoidal  layers  (5+10%,L=5,H=2.5) 

10 

0.083 

-0.007 

0.090 

2-Gaussian  layer  (4.5+10%/5.0+10%,2km) 

11 

-0.008 

-0.048 

0.040 

steeply  dipping  layers  (52th) 

12 

-0.024 

-0.057 

0.033 

gently  dipping  layers  (26°) 

13 

0.015 

-0.086 

0.101 

steeply  dipping  layers  (-52°) 

14 

-0.001 

-0.103 

0.102 

gently  dipping  layers  (-26°) 
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Figure  7.  P  wave  in  a  half  space  (a  =  6.0  km/s,  p  =  3.SS  km/s)  incident  at  20s  upon  a  2  km  layer  with 
average  P-wave  velocity  of  5  km /s  and  a  self-similar  10%  rms  velocity  variation  superimposed  by  a 
gentle  topography  (indicated  in  the  0  sec  frame,  also  model  2  in  Tables  1  and  2)  The  S-wave  velocity  is 
assumed  to  be  proportional  to  the  P-wave  velocity.  Darkness  of  the  snapshots  are  proportional  to  the 
displacement  amplitude.  Snapshots  of  the  displacement  field  are  shown  at  1  second  intervals.  The  dila- 
tational  strain  is  recorded  at  32  locations  at  a  depth  of  O  S  km  in  order  to  infer  the  excitation  of  far-field 
P  waves  from  explosion  sources.  Although  absorbing  boundary  conditions  are  used,  care  must  be  taken 
to  avoid  residual  reflections  from  the  sides  of  the  grid. 
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Figure  9.  Synthetic  far-field  P-  (top)  and  SV-wave  (bottom)  inferred  by  the  principle  of  reciprocity  for 
model  2.  The  original  dilatational  strain  history  (5  Hr  low-pass)  responding  to  incident  broadband  P  or 
SV  plane  wave  recorded  at  32  locations  at  0.5  km  depth  in  the  reference  model.  Shown  here  are  the 
deconvolved  synthetics  corresponding  to  VSB  50  KT  in  hard  rock.  The  peak  amplitude  of  these  syn¬ 
thetics  was  measured  and  compared  to  the  average  peak  amplitude  of  the  reference  model. 
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Figure  10.  Average  spectral  ratio  as  a  function  of  frequency.  Log(P\IPo)  (upper)  and  Log(.Lg]lLg0 ) 
(lower),  of  the  Model  2  relative  to  the  reference  model.  P  wave  response  of  model  l  in  the  0.5  to  1 .0 
Hz  range  is  deficient  with  respect  to  the  reference  model  by  0.348  log  units,  while  the  S  wave  response 
is  0.063  log  unit  above  the  reference  model.  Vertical  bars  represent  the  standard  error  of  a  single 
observation. 
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Figure  11.  Averaged  P  and  SV(Lg)  peak  amplitudes  of  array  of  shallow  explosions  in  fourteen  crustal 
models  sorted  with  the  \ogi0(Lg/Lg0)  values.  Several  observations  are  obvious:  (1)  Dipping  layers 
(models  1 1  through  14)  generate  smaller  Lg  than  the  normalizing  model,  while  they  all  generate  more  P 
than  the  reference  model.  (2)  Media  with  topography  (e.g.  models  1  through  4)  which  represent 
CEKTS  all  generate  more  Lg  than  the  normalizing  model,  while  they  excite  less  P  due  to  strong  P  to  S 
conversion.  (3)  Dipping  layers  (models  11  through  14)  are  more  efficient  than  ail  other  models  for  P 
excitation.  Thus  m*(P)  and  m»(Lg)  appear  to  be  negatively  correlated. 
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Figure  12.  Same  as  Figure  1 1  except  the  P  and  SV(Lg)  excitations  are  measured  with  the  averaged 
spectra!  content  in  [0.3,1]  Hz  band.  Crustal  models  with  topography  generate  more  Lg  and  less  P  than 
models  with  dipping  layers,  same  as  the  result  derived  from  peak  amplitude  measurement. 
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E.  Kazakh  Events,  CDSN-WMQ 

Figure  13.  Short  period  seismograms  of  two  Shagan  events  87171  (78.74E,  49.91N,  mb-6.1)  and 
87347  (78.8SE,  49.96N,  mb-6.1).  and  a  Degelen  event  87198  (78.1  IE  49.80N,  mb-5.8)  recorded  at 
CDSN  station  WMQ.  Each  trace  is  scaled  by  the  peak  amplitude.  Note  the  relatively  less  P  energy 
(with  respect  to  Lg  energy)  in  the  Degelen  event  87198  as  compared  to  Shagan  events  of  similar  mag¬ 
nitudes.  This  observation  is  consistent  with  the  finite-difference  experiments  (c/..  Tables  1,  2  and  Fig¬ 
ures  11,  12). 
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DISCUSSION  AND  SUGGESTIONS 

We  see  that  moderate  heterogeneity  in  a  half  space  does  not  attenuate  short-period 
fundamental  Rayleigh  waves  nearly  so  much  as  rough  topography  does,  but  it  can  still 
contribute  substantial  P-coda  and  moderate  attenuation  of  Rg.  For  the  Gaussian  media 
used  in  these  simulations,  the  energy  lost  due  to  body  wave  conversion  varies  from 
several  percent  to  20%  in  12  km  distance.  A  significant  result  of  the  simulations  is  that 
reflection  of  Rayleigh  waves  by  heterogeneity  at  normal  incidence  is  in  most  cases 
inefficient,  as  was  the  case  for  rough  topography.  The  only  exception  observed  was  a 
folded  structure  with  a  resonant  response  to  the  incident  Rayleigh  wave.  Therefore  we 
should  not  expect  to  see  Rayleigh-wave  back-scattering  as  a  significant  contributor  to  the 
multiple  scattering  of  fundamental  Rayleigh  waves  that  populate  coda  for  Gaussian  or 
self-similar  media.  Backscattering  can  be  significant  for  media  with  well  defined  organ¬ 
ized  structures  such  as  folded  sedimentary  structures.  In  such  a  case  the  backscattered 
wave  has  a  narrow  narrow  bandwidth  reflecting  the  resonance  phenomenon. 

More  complicated  random  media  contain  many  scale  lengths  and  introduce  broad¬ 
band  effects.  The  Rayleigh-wave  attenuation  is  a  complicated  function  of  frequency,  but 
it  attains  a  maximum  in  the  range  where  the  characteristic  wavelength  of  the  medium 
matches  the  wavelength  of  the  incident  Rayleigh  wave. 

At  low  frequency  (X  >  a)  the  coda  dilatational  and  rotational  wavenumber  spectra 
indicate  that  the  scattered  P  and  SV  waves  are  scattered  in  the  forward  direction  except 
for  strongly  folded  structures  or  Gaussian  medium  with  very  strong  velocity  fluctuations 
(20%).  For  higher  frequencies  (X  <  a),  the  scattered  body  waves  are  always  maximum  in 
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the  forward  direction.  A  detailed  analysis  of  Rayleigh  to  P-coda  scattering  will  have  to 
take  this  effect  into  account  with  an  effective  radiation  pattern  to  the  "equivalent 
scatterer". 

The  results  presented  in  AFGL-TR-87-0322  offer  a  beginning  approach  to  under¬ 
standing  the  effects  upon  surface  waves  of  scattering  by  lateral  heterogeneity.  A  com¬ 
plete  exploration  of  the  problem  will  require  variation  of  the  P  and  S  velocities,  near  sur¬ 
face  velocity  gradients,  and  crustal  velocity  heterogeneity  as  a  function  of  depth.  Much  of 
the  SV  energy  scattered  by  near  surface  heterogeneities  is  concentrated  at  apparent  veloci¬ 
ties  within  150%  of  the  Rayleigh  phase  velocity.  A  crust  with  surface  layer  with  P  =  2.96 
km/sec,  as  was  used  in  these  simulations,  would  leak  much  of  the  energy  to  the  mantle. 
However,  if  the  near-surface  velocity  of  the  model  is  lowered,  then  the  slowness  space 
occupied  by  the  scattered  waves  will  scale  to  the  Rayleigh  phase  velocity,  and  more 
energy  will  be  trapped  in  the  crust.  Other  ways  of  increasing  the  trapping  of  the  scattered 
SV  energy  are  increasing  the  oc/p  ratio  in  the  near  surface,  introduction  of  gradients  near 
the  surface  to  create  higher  order  modes,  and  introduction  of  deeper  velocity  hetero¬ 
geneity  to  scatter  P  and  SV  energy  back  into  the  waveguide.  In  short,  all  these  mechan¬ 
isms  can  act  only  to  increase  the  Rayleigh-to-Lg  coupling.  Therefore,  we  expect  that  in 
real  seismological  situations  much  of  the  Rayleigh-wave  energy  scattered  into  SV  by 
near-surface  heterogeneity  will  be  trapped  in  the  crust  and  will  find  a  path  to  the  Lg 
wavepacket. 

In  our  simplistic  attempt  to  model  SV-Lg  excitation  due  to  near-source  heterogeneity 
scattering  effects  (AFGL-TR-88-0025 ),  we  found  that  the  total  variance  of  the  spectra  of 
the  teleseismic  P  wave  was  greater  than  the  variance  of  the  P-SV-Lg  excitation  (with  the 
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slowness  we  investigated).  Although  we  found  that  different  geologic  models  gave 
different  coupling,  the  variance  was  always  larger  for  P  than  for  P-SV  for  any  models. 
This  may  explain  Nuttli’s  claim  that  Lg  is  a  stable  estimator  in  a  fixed  geologic  setting. 

While  we  continue  to  experiment  with  various  models,  our  preliminary  results  indi¬ 
cate  that  P  to  SV  conversion  is  strongly  enhanced  by  velocity  variation  in  the  vicinity  of 
rough  topography  and  by  the  introduction  of  low  velocity  layers  near  the  surface.  The 
introduction  of  interfaces  alone  does  not  of  itself  increase  SV  excitation  with  the  required 
slowness  range.  We  continue  to  experiment  with  the  geometry  of  heterogeneity  and  with 
the  scale  lengths  of  the  heterogeneity.  Although  we  cannot  presently  explain  Nuttli’s 
(1987)  results,  we  predict  substantial  variations  in  SV  Lg  excitation  by  explosions  embed¬ 
ded  in  crustal  heterogeneity.  It  seems  that  P-to-SV  is  not  the  only  mechanism  for  explo¬ 
sion  Lg  excitation,  so  it  is  necessary  to  investigate  the  excitation  of  SH(Lg)  as  well. 

A  possibility  is  that  Nuttli’s  mb:Lg  relationship  might  be  related  to  Rayleigh-to-P 
conversion  away  from  the  immediate  location  of  the  source.  Our  numerical  simulations 
treated  only  the  P-S  conversions  that  might  occur  within  a  few  km  of  the  source,  given 
some  simple  models.  If  either  the  Rayleigh  excitation  or  the  Rayleigh  scattering  is 
different  for  CEKTS  and  EEKTS,  then  we  could  see  the  difference  in  Rayleigh  to  Lg. 
Since  the  two  locations  are  only  70  km  apart,  the  Rayleigh-to-Lg  difference  would 
presumably  have  to  occur  in  the  first  20-25  seconds.  Thus  it  seems  necessary  to  examine 
whether  the  P-coda  are  different  for  Degelen  and  Shagan  in  the  first  20-25  seconds. 
Similarly,  P-SV  conversion  could  be  happening  further  away  from  the  source  than  we  are 
modeling.  It  is  also  possible  that  the  non-linear  source  effects  might  produce  larger  SV  at 
one  site  versus  another.  These  hypotheses  as  well  as  the  3-dimensional  effects  were  not 
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addressed  in  our  current  experiments. 

As  recommendations  for  further  work  along  these  lines,  we  would  suggest  the  fol¬ 
lowing  tasks: 

(1)  Extensions  of  the  current  finite-difference  code  from  2-D  to  3-D  to  study  the 
attenuation  of  body  waves  by  3-D  heterogeneity  in  the  crust,  to  test  hypotheses 
about  the  generation  of  P  coda  and  anisotropic  P  wave  generation,  and  to  model  the 
generation  of  transverse  Lg  by  explosions. 

(2)  Introduction  of  other  numerical  methods  to  explore  the  coupling  (scattering)  of 
modes  of  wave-guide  regional  phases  such  as  Pg  and  Lg,  as  well  as  die  scattering  of 
Pn  and  Sn.  These  methods  include  2-D  and  3-D  scattering  from  localized  hetero¬ 
geneity  as  well  as  from  rough  boundaries. 

(3)  Coupling  of  efficient  reflectivity  methods  to  finite  difference  calculations  to  pro¬ 
pagate  the  scattered  field  to  regional  distances  and  to  drive  the  finite  difference 
responses  with  realistic  in-coming  regional  phases. 

(4)  Investigation  of  scattering  of  fundamental  and  higher  mode  short-period  Rayleigh 
waves  by  2-D  topography  and  shallow  heterogeneity  with  more  realistic  velocity  gra¬ 
dients  near  the  surface. 

(5)  Extension  of  the  general  topographic  boundary  condition  to  include  the  general 
fluid-solid  interface  for  the  modeling  of  scattering  at  rough  fluid-solid  boundaries. 

(6)  Improvement  of  the  polygonal  free-surface  boundary  conditions  for  higher  precision. 
One  distinguishable  feature  of  our  FD  code  is  that  it  allows  simulations  with  fairly 
rough  topography.  The  algorithm  we  use  (Jih  et  al.,  1988)  has  first-order  accuracy 
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consistent  with  the  standard  one-sided  extrapolation  formula  for  the  flat  free-surface. 
Even  though  a  number  of  FD  techniques  have  been  proposed  in  recent  years  such  as 
using  higher-order  spatial  difference  operator  at  the  interior  points,  implicit  rather 
than  explicit  scheme,  etc.,  none  of  them  have  demonstrated  significant  improvement 
over  the  traditional  explicit  second-order  scheme  as  long  as  the  only  available 
scheme  implementing  irregular  free-surface  has  accuracy  only  of  order  one.  This  is 
because  the  overall  accuracy  of  a  FD  scheme  would  degrade  if  the  boundary  condi¬ 
tions  are  represented  by  a  scheme  with  accuracy  lower  than  that  for  the  interior 
medium.  We  suggest  deriving  an  improved  FD  formulation  of  the  boundary  condi¬ 
tions  so  that  the  accuracy  can  be  compatible  with  at  least  second  order  in  spatial 
increment. 
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Section  A.O 
SUMMARY 

This  first  volume  of  the  user's  manual  gives  a  detailed  description  of  the  current 
version  of  a  sofware  package,  TGALFD,  for  generating  synthetic  seismograms.  The 
program  has  been  developed  at  Teledyne  Geotech  Alexandria  Laboratories  (TGAL) 
during  the  past  several  years.  This  package  consists  mainly  of  a  2-dimensional  2nd- 
order  explicit  finite-difference  (FD)  code  which  permits  various  source  types,  topo¬ 
graphical  free  surface,  as  well  as  arbitrary  fluctuation  in  (2-dimensional)  medium  pro¬ 
perties.  Some  sample  runs,  free- surface  boundary  conditions  used  for  topography  han¬ 
dling,  and  a  listing  of  the  source  code  are  included.  The  supporting  routines  used  in 
analyzing  the  output  synthetics,  routines  used  for  media  generation,  as  well  as  the 
(machine-dependent)  color  snapshots  display  routines  will  be  described  in  a  follow  up 
volume  of  this  manual. 

This  FORTRAN-77  code  has  been  run  under  the  UNIX  operation  system  on 
VAX,  SUN,  Celerity,  and  Convex  computers.  Users  of  the  Center  for  Seismic  Studies 
may  contact  Geotech  for  the  use  of  this  package. 
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Section  A.l 

FINITE  DIFFERENCE  METHOD 

Wave  propagation  problems  in  seismology  involve  the  solution  of  a  set  of 
differential  equations  in  a  medium  in  which  the  material  properties  vary  in  space,  i.e., 
in  the  earth.  The  use  of  numerical  simulations  is  a  straightforward  means  for  studying 
this  kind  of  problem,  especially  when  laterally-varying  velocity  structure  or  complex 
topographic  relief  is  encountered.  Methods  such  as  Gaussian  beam  technique  and  ray 
theoretical  schemes  are  restricted  to  cases  where  variations  of  the  medium  are  much 
larger  than  the  seismic  wavelength.  The  Kirchhoff-Helmholtz  integration  method  is 
useful  for  media  with  sharp  interfaces,  but  it  doesn’t  include  the  multiple  scattering 
along  the  interfaces,  and  it  is  not  appropriate  for  reflections  from  velocity  gradients 
similar  in  extent  to  the  seismic  wavelength.  Perturbation  methods  arc  applicable  only 
for  weak  scattering  problems.  Among  all  numerical  approaches,  finite-difference  (FD) 
and  finite-element  (FE)  methods  are  not  restricted  to  velocity  variations  of  a  particular 
size  with  respect  to  wavelength.  FD  and  FE  can  generate  synthetic  seismograms  for 
very  complicated  media  in  cases  of  weak/strong  or  multiple  scattering. 

FD  method  solves  either  the  wave  equations  or  the  elastodynamic  equations  by 
replacing  the  partial  derivatives  in  space  and  time  by  their  FD  approximations.  When 
explicit  FD  method  is  used,  which  is  the  most  popular  FD  technique  to  date,  the 
wavefield  of  a  specific  time  instant  is  solved  one  grid  point  by  one  grid  point  with 
nearby  grid  data  at  previous  steps.  For  schemes  that  use  second-order  approximations 
to  the  temporal  derivative,  only  two  grid  planes  of  displacement  (or  stress,  velocity) 
must  be  stored  to  perform  the  updating  process.  Once  the  entire  grid  is  updated,  FD 
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then  proceeds  to  compute  the  wavefield  of  next  time  instant  until  a  certain  preselected 
number  of  iterations  is  reached.  The  output  of  FD  method  can  be  snapshots  of  the 
entire  grid  at  specific  times  or  synthetic  seismograms  recorded  at  specific  grid  points. 

The  excellent  review  papers  by  Chin  et  al.  (1984),  Frankel  (1988),  and  Stephen 
and  Burton  (1988)  contain  more  detailed  discussion  of  finite-difference  method  as  well 
as  other  numerical  methods  used  in  seismology.  A  fairly  comprehensive  list  of  refer¬ 
ences  is  given  at  the  end  of  this  manual. 
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Section  A.2 

TGAL’S  PROGRAM  FD8 

Revision  History 

Teledyne  Geotech  has  been  engaged  in  the  development  and  utilization  of  FD 
code  for  a  long  time.  Z.  A.  Der,  J.  Bumetti,  and  T.  McElfresh  initialized  the  code 
design  in  which  tne  2nd  order  explicit  FD  formulation  (Kelly  et  al.,  1976)  was 
adopted.  They  used  a  monochromatic  P/SV  planar  source  as  well  as  symmetric  boun¬ 
dary  conditions  to  model  the  effects  of  crustal  structure  on  teleseismic  and  some 
regional  phases  during  1978-1981  (Der  etal.,  1978;  Barker  etal.,  1981).  K.  L. 
McLaughlin,  T.  McElfresh,  and  L.  Anderson  implemented  an  Ohnaka  (broadband) 
P/SV  source,  a  point  (line)  source,  and  absorbing  boundary  conditions  (Clayton  and 
Engquist,  1977;  Emerman  and  Stephen,  1983)  during  1983-1985.  R.-S.  Jih  developed 
the  lst-order  representation  of  the  free-surface  boundary  condition  to  handle  polygonal 
topography  (Jih  et  al.,  1988)  and  coded  a  source-independent  fundamental  Rayleigh 
wave  generation  routine  (Boore,  1970;  Munasinghe  and  Famell,  1973;  Levander, 
1985). 

The  current  version  of  FD  code  fd8  is  quite  different  from  all  earlier  versions 
(fdabci  through  fdabc6),  after  a  series  of  heavy  revisions  was  performed  during 
1986-1987  to  allow  more  options  to  model  realistic  problems,  even  though  several  sub¬ 
routines  still  retain  their  original  names.  This  code  has  been  utilized  extensively  by 
researchers  at  Alexandria  Laboratories  to  study  various  seismological  problems  with 
funding  from  DARPA.  The  following  discussion  will  therefore  be  confined  to  fd8 
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fd8  reads  in  control  parameters  from  the  standard  input,  and  seismograms  (time 
series)  and  error  messages  are  written  to  the  standard  output.  Snapshots  of  the 
horizontal/vertical  displacements  and/or  strains  are  stored  in  the  output  file  specified  by 
the  input  file. 


Initial  Conditions 


The  initial  wave  could  be 


(+1)  broadband  planar  P  wave  of  Ohnaka  shape, 

(-1)  broadband  planar  S  wave  of  Ohnaka  shape, 

(+2)  monochromatic  planar  P  wave  of  sinusoidal  shape, 

(-2)  monochromatic  planar  S  wave  of  sinusoidal  shape, 

(+3)  pure  compressional  (P)  wave  generated  at  a  single  point, 

(-3)  double  couple  (S)  point  source, 

(+4)  fundamental  mode  Rayleigh  wave  with  Ricker  wavelet  shape, 

(+5)  a  (time  series)  driver  file  shaking  a  single  point, 

(+6)  arbitrary  wavefield  setting, 

(+7)  broadband  planar  P  wave  of  Gaussian  shape, 

(-7)  broadband  planar  S  wave  of  Gaussian  shape, 

Except  for  options  (+5)  and  (+6)  in  which  the  source  file  or  initial  wavefields 
must  be  generated  elsewhere  in  advance,  fd8  is  completely  self-contained  to  initialize 


the  wavefields  at  2  consecutive  time  instants  to  generate  proper  wave  propagation  later 
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Boundary  Conditions 

Absorbing  boundary  condition  is  the  default  for  the  bottom  and  side  edges.  Sym¬ 
metric  side  boundary  conditions  are  used  only  for  the  case  of  normally  incident  planar 
waves.  Free  surface  is  assumed  on  the  top  of  the  grid  whenever  a  nontrivial  topogra¬ 
phy  is  involved.  All  these  boundary  conditions  can  be  altered  by  choosing  appropriate 
input  control  parameters.  For  instance,  in  the  case  of  a  point  source  (i.e.  option  +/-3 
or  5)  without  topography,  there  are  3  more  choices  by  playing  with  incident  angle: 

(1)  O-degree  causes  all  4  edges  to  be  absorbing, 

(2)  360-degree  causes  symmetric  top  plus  absorbing  sides,  bottom. 

(3)  720-degree  causes  all  4  edges  to  be  symmetric. 

These  extra  options  are  meant  mainly  to  demonstrate  the  effects  of  miscellaneous 
boundary  conditions  rather  than  to  model  realistic  seismological  problems. 

One  distinguishable  feature  of  fd8  is  that  it  allows  simulations  with  fairly  rough 
topography.  The  algorithm  (Jih  et  al. ,  1988)  used  here  is  an  improved  version  of  Ilan 
(1977).  On  the  inclined  firee-surface,  the  vanishing  stress  conditions  are  implemented  to 
a  rotated  coordinate  system  parallel  to  the  inclined  boundary,  as  previous  works  did. 
For  each  transition  point  on  the  topography  where  the  slope  changes,  we  use  the  first- 
order  approximation  of  boundary  conditions  in  a  locally  rotated  coordinate  system  in 
which  the  normal  axis  always  coincides  with  the  bisector  of  the  comer.  These  extra¬ 
polation  formulae  are  consistent  with  boundary  conditions  to  first-order  accuracy  in 
spatial  increment,  same  as  the  classical  one-sided  explicit  approximation  scheme 
widely  used  for  flat  free-surface  case.  Testing  results  indicate  that  this  scheme  works 
stably  for  fairly  complicated  geometric  shapes  consisting  of  ridges  and  valleys  with 
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steep  and  gentle  slopes  over  a  range  of  Poisson  ratios  of  practical  interest,  thus  ena¬ 
bling  us  to  study  more  realistic  problems  for  which  the  topography  plays  a  significant 
role  in  shapu.^  he  wavefield  and  for  which  an  analytical  solution  might  not  be  avail¬ 
able. 

Output 

The  program  converts  numerical  wavefields  into  character  wavefields  and  stores 
these  snapshots  in  an  ASCII  text  file.  The  output  wavefield  could  be  the  whole  grid  or 
only  the  central  portion.  An  input  parameter  determines  whether  fixed  gain  or 
automatically  adjustable  scale  is  used  in  the  conversion  procedure. 

Displacements  and/or  strain  may  be  recorded  as  time  series  at  any  interior  points 
for  the  strain  or  any  grid  points  for  the  displacements. 

The  program  also  stores  the  wavefields  at  2  consecutive  instants  and  all  required 
parameters  at  a  prespecified  rate  so  that  it  can  be  re-started  in  case  the  job  is  ter¬ 
minated  in  the  middle. 
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Section  A.3 

SAMPLE  INPUT  FILES  AND  SAMPLE  RUNS 


Sample  Input  Parameter  File  for  fd8 

grid  dimension:  kw,kh 
100  100 

x,z  spacings  of  the  grid  mesh  &  temporal  spacing:  dx,dz,dt 

2  2  .2 

water  level:  iwater  (must  be  >  2) 

500 

homogeneous  or  heterogeneous  medium  flag  (0  or  1) 

1 

inLm 

inMu 

inRho 

topography  (sea-floor  or  ground)  model:  inTOP 
TOP 

file  name  for  the  output  snapshots 
movie 

choose  wave  type:  itype  (see  Section  A.2  for  legal  options) 

7 

incidence  angle  (degree),  iOjO,  wavelength  (km)  (see  Remarks  (a),  (b)) 

0.01  71  2  20 

option  for  snapshot  display:  component  flag,  AGC  flag,  whole  (Remark  (c)-(e)) 
#  of  stations  to  record  seismograms 
10 

coordinates  of  sensors:  (Remarks  (0) 

071  171  271  371  471  571  671  771  871  971 

02  02  02  02  02  02  02  02  02  02 

total  number  of  iterations;  iterations  per  snapshot  generation 

3000  200 


Remarks: 

(a)  Wavelength  is  needed  for  point  source,  sine  or  Rayleigh  wave, 

(b)  Location  of  source  can  also  be  used  to  adjust  boundary  conditions, 

(c)  6  options  for  snapshot  displaying  are  available 

h  :  horizontal  wave  fields, 
v  :  vertical  wave  fields, 
b  :  both  horizontal  &  vertical  wave  fields, 
s  :  dilatational  &  rotational  strain  fields, 
m  :  combined  displacement  fields. 
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(d)  AOC  >  0  =>  resolution  adjusted  for  each  frame  individually, 

(e)  size  flag  >  0  =>  whole  grid  will  be  shown, 

(f)  List  X-coordinates  of  all  sensors  first  Negative  X-coordinate  means  strain  sensor. 
Y-coordinates  should  be  in  another  record,  and  should  be  consistent  with  #  of 
sensors  above. 


Sample  Topography  File 
10 

00  00  00  01  02  03  04  05  06  07  08  08  08  08  07  06  05  05  05  05 
EOF 

Remarks: 

(g)  The  1st  line  gives  the  reference  floor  level  of  the  topography  (counted  from  the 
top  of  the  grid).  The  2nd  line  gives  elevation  with  respect  to  the  reference  level 
at  each  column,  e.g. ,  03  means  free  surface  is  at  7th  (=  10  -  3)  row  in  Z-direction 
(counted  from  above).  Note  that  the  reference  level  must  be  deep  enough,  and 
each  segment  of  the  polygon  must  contain  at  least  2  sub-segments  before  the 
slope  changes. 


Sample  Density  Model  File 
grid  size:  20  25 

grid  spacing:  0.1000  0.1000  km 

Self-Similar  Medium 
Extracted  from  another  model 

5  (dummy  line) 

6  (dummy  line) 

7  (dummy  line) 

8  (dummy  line) 

0.4927E+01  0.4982E+01  0.5106E+01  0.5123E+01  0.5128E+01 

0.5078E+01  0.5213E+01  0.5393E+01  0.5229E+01  0.5239E+01 

0.5584E+01  0.5685E+01  0.5583E+01  0.5542E+01  0.5256E+01 


EOF 

Remarks: 

(h)  Lines  1  and  2  specify  the  grid  size  and  spacings,  and  lines  3  and  4  are  for 
identification  purpose.  Line  5  thru  8  are  dummy.  The  remaining  lines  give 
Lame’s  constant  at  (i,j),i=l,kw),j=l,kh,  X,  p.  fields  have  the  same  format  as  p 
field. 
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Just  as  a  simple  demonstration  of  the  capabilities  as  well  as  the  limitations  of  our 
code,  figures  1  through  4  give  the  snapshots  of  wave  propagation  with  various  optional 
conditions  imposed. 

Example  A.1 

Figure  A.l  shows  the  propagation  of  a  normally  incident  plane  P  wave  through  a 
model  with  a  45°  ramp  on  the  top  of  grid  and  von  Neumann  (i.e.  0-slope  or  sym¬ 
metric)  boundary  condition  used  on  both  sides.  The  appropriate  P-S  conversions  and 
the  reflections,  the  diffractions  satisfying  Snell’s  law  and  Huygen’s  principle  are 
clearly  visible  in  these  successive  snapshots  taken  every  second.  It  is  easy  to  verify 
the  first-order  accuracy  in  spatial  increment  of  our  one-sided  explicit  representation  of 
the  free-surface  boundary  conditions  in  this  case. 

Example  A.2 

Figure  A.2  shows  the  snapshots  of  displacement  fields  generated  by  an  upgoing  P 
wave  in  a  grid  with  steep  topographic  configuration.  The  successive  frames  separated 
by  0.125  sec  show  the  initialization  of  the  wave  (A),  P-reflection  followed  by  S  wave 
starting  at  right  (B,C),  completely  developed  reflections  from  all  parts  of  the  topogra¬ 
phy  (E)  «nd  complex  wavefields  containing  reflections,  diffractions,  and  possibly 
excited  surface  waves  (E,F).  The  initial  P  wave  has  an  incident  angle  of  20°  and  the 
topography  is  a  (due  north  344°)  cross-section  of  Taourirt  Tan  Afella  Massif  in  south¬ 
ern  Algeria.  It  can  be  observed  that  the  free-surface  reflection  is  severely  altered  due 
to  scattering  from  the  free-surface.  The  symmetric  boundary  conditions  on  the  sides 
cause  some  undesirable  edge  reflections  from  the  left  side  at  a  later  stage  of  the 


Final  Report 


37 


March  1988 


AFGL-TR-W-0093 


TGAL4S42 


simulation. 

Example  A.3 

Figure  A.3  uses  the  same  topographical  configuration  as  in  Figure  A.2  with 
absorbing  boundary  conditions  (Clayton  and  Engquist,  1977;  Emerman  and  Stephen, 
1983)  adopted  on  the  sides  and  bottom  to  suppress  the  artificial  reflections  from  the 
sides  of  the  grid.  The  compressional  point  source  in  this  2-D  rectangular  scheme  is  in 
fact  a  line  source  in  3-D  space  and  hence  is  not  realistic  enough  in  some  aspects. 
Note  that  the  quasi-transparent  boundary  conditions  allow  the  wave  to  disappear  into 
the  sides  and  bottom  of  the  grid. 

Example  A.4 

Figure  A.4  shows  a  Rayleigh  wave  incident  on  a  rough  topographic  profile  super¬ 
imposed  on  a  grid  with  absorbing  boundary  conditions  for  the  sides  and  the  bottom. 
Figures  (A)  through  (F)  correspond  to  displacement  wavefields  at  distinct  times  with  a 
temporal  spacing  of  2  sec.  Note  that  the  high  frequency  scattering  of  Rayleigh  wave 
is  forward  (McLaughlin  and  Jih;  1986). 
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Section  A.4 
REFERENCES 

For  the  user’s  convenience,  references  are  hereby  divided  into  three  categories: 
(1)  publications  directly  used  in  coding  fd8,  (2)  TGAL’s  research  that  utilized  fd8  or 
its  earlier  version,  and  (3)  general  references.  All  Teledyne  Geotech  reports  are  avail¬ 
able  through  the  National  Technical  Information  Service. 
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Figure  A.1  The  propagation  of  a  normally  incident  plane  P  wave  through  a  model  with  a  45"  ramp  on 
the  top  of  grid  and  symmetric  boundary  condition  used  on  both  sides.  The  appropriate  P-S  conversions 
and  the  reflections,  the  diffractions  satisfying  Snell’s  law  and  Huygen’s  principle  are  clearly  visible  in 
these  successive  snapshots  taken  every  second  (from  Jih  et  ai,  1988). 
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Figure  A.2  The  displacement  fields  generated  by  a  plane  P  wave  of  incidence  angle  20s  in  a  grid  with 
steep  topographic  configuration.  The  successive  frames  separated  by  0.12S  sec  show  the  initialization 
of  the  wave  (A),  P-reflecdon  followed  by  S  wave  starting  at  right  (B,C),  completely  developed 
reflections  from  all  parts  of  the  topography  (E)  and  complex  wavefields  containing  reflections, 
diffractions  and  possibly  excited  surface  waves  (E,F).  It  can  be  observed  that  the  free-surface  reflection 
is  severely  altered  due  to  scattering  from  the  free-surface  (from  Jih  et  al„  1988). 
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(E)  (F) 

Figure  A .3  Same  topographical  configuration  as  in  Figure  A.2  with  a  compressional  point  (line)  source, 
and  absorbing  boundary  conditions.  Note  that  the  quasi-transparent  boundary  conditions  allow  the  wave 
to  disappear  into  the  sides  and  bottom  of  the  grid.  Snapshots  are  separated  by  0.23  second. 


Final  Report 


49 


March  1988 


APGL-m*W093 


TGAL4M-C2 


(E)  (F) 

Figure  A.4  Rayleigh  wave  incident  on  a  rough  topographic  profile  superimposed  on  a  grid  with  absorb¬ 
ing  boundary  conditions  for  the  sides  and  the  bottom.  Figures  (A)  through  (F)  correspond  to  displace¬ 
ment  wavefields  at  distinct  times  with  a  temporal  spacing  of  2  sec.  Note  that  the  high  frequency 
scattering  of  Rayleigh  wave  is  forward  (McLaughlin  and  Jih;  1986). 
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